// panda.cpp : Defines the entry point for the console application.
//

#if _MSC_VER > 1000
#pragma warning( disable :  4511 4512 4786)
// nasty hack for VC++ non-compliance with standard C++ for loop declaration scope
#define for if (0); else for
#endif

#include <stdio.h>
#include <math.h>
#include <string>
#include <set>
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <iterator>
#include <sstream>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <boost/tuple/tuple.hpp>
#define _MY_STL_MACROS
#include "components/stlhelper.h"
#include "components/perlutil.h"
#include "components/progress_indicator.h"
#include "components/ultoa.h"
#include "components/pg_copy_guard.h"
#include "components/pg_db_cursor.h"
#include "components/portable_timer.h"
#include <assert.h>
#include "get_arg.h"

using std::string;
using std::sort;
using std::stringstream;
using std::cerr;
using std::cin;
using std::vector;
using std::fill;
using std::cout;
using std::ifstream;
using std::ofstream;
using std::ostream;








//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//	loci code
//
struct t_segment
{
	unsigned beg;
	unsigned end;
	t_segment(unsigned beg_, unsigned end_)
		: beg(beg_), end(end_){}
	t_segment(const string& data)
	{
		beg = atoi(data.c_str()) - 1;
		end = atoi(data.c_str() + data.find(':') + 1) - 1;
		if (beg  > end)
			beg = end = 0;
	}
};
typedef vector<t_segment> t_loci;



//________________________________________________________________________________________

//	merge_overlapping_loci

//________________________________________________________________________________________
void merge_overlapping_loci(t_loci& loci,
							t_loci& complement_loci,
							unsigned aln_beg, unsigned aln_end)
{
	vector<char> bits(aln_end, '\x0');
	typedef vector<char>::iterator iter_pos;

	// mark bits of all masked regions
	for (unsigned i = 0; i < loci.size(); ++i)
	{
		fill(	bits.begin() + std::min(loci[i].beg, aln_end),
				bits.begin() + std::min(loci[i].end + 1, aln_end),
				0x1);
	}

	// non-overlapping masked loci
	t_loci loci2;

	// get masked regions
	iter_pos first = find(bits.begin() + aln_beg, bits.end(), 0x1);
	while (first != bits.end())
	{
		iter_pos last = find(first+1, bits.end(), 0x0);
		loci2.push_back(t_segment(first - bits.begin(),
								   	   last  - bits.begin() - 1));
		first = find(last, bits.end(), 0x1);
	}
	loci2.swap(loci);

	// get non-masked regions
	first = find(bits.begin() + aln_beg, bits.end(), 0x0);
	while (first != bits.end())
	{
		iter_pos last = find(first+1, bits.end(), 0x1);
		complement_loci.push_back(t_segment(first - bits.begin(),
												 last  - bits.begin() - 1));
		first = find(last, bits.end(), 0x0);
	}
}

//________________________________________________________________________________________

//	subtract_overlapping_loci

//________________________________________________________________________________________
void subtract_loci(	const t_loci& original_loci, const t_loci& remove_loci,
					t_loci& resulting_loci,
			   		unsigned aln_beg, unsigned aln_end)
{
	resulting_loci.clear();
	vector<char> bits(aln_end + 1, '\x0');
	typedef vector<char>::iterator iter_pos;

	// mark bits of all loci
	for (unsigned i = 0; i < original_loci.size(); ++i)
	{
		fill(	bits.begin() + std::min(original_loci[i].beg, aln_end),
				bits.begin() + std::min(original_loci[i].end + 1, aln_end),
				0x1);
	}

	// remove bits of all loci
	for (unsigned i = 0; i < remove_loci.size(); ++i)
	{
		fill(	bits.begin() + std::min(remove_loci[i].beg, aln_end),
				bits.begin() + std::min(remove_loci[i].end + 1, aln_end),
				0x0);
	}

	// get resulting loci
	iter_pos first = find(bits.begin() + aln_beg, bits.end(), 0x1);
	while (first != bits.end())
	{
		iter_pos last = find(first+1, bits.end(), 0x0);
		resulting_loci.push_back(t_segment(first - bits.begin(),
			   						   last  - bits.begin() - 1));
		first = find(last, bits.end(), 0x1);
	}
}


//________________________________________________________________________________________

//	 extract_seq_pair_with_loci

//________________________________________________________________________________________
inline string extract_seq_with_loci(const std::string& seq, const t_loci& loci)
{
	string nrvo;
	for (unsigned i=0; i < loci.size(); ++i)
		nrvo += seq.substr(loci[i].beg,  loci[i].end + 1 - loci[i].beg);
	return nrvo;
}
inline void extract_seq_pair_with_loci(vector<string>& sequences,
										const std::string& seq1,
										const std::string& seq2,
										const t_loci& loci)
{
	// add sequence for mouse
	string seq_region;
	for (unsigned i=0; i < loci.size(); ++i)
		seq_region += seq1.substr(loci[i].beg,  loci[i].end + 1 - loci[i].beg);
	sequences.push_back(seq_region);

	// add sequence for human
	seq_region.clear();
	for (unsigned i=0; i < loci.size(); ++i)
		seq_region += seq2.substr(loci[i].beg,  loci[i].end + 1 - loci[i].beg);
	sequences.push_back(seq_region);
}

//________________________________________________________________________________________

//	 extract_cdna_pair_with_loci

//________________________________________________________________________________________
struct op_prot_loci_to_dna
{
	void operator() (t_segment& a)
	{
		a.beg *= 3;
		a.end = (a.end + 1) * 3 - 1;
	}
};

inline void extract_cdna_pair_with_loci(vector<string>& sequences,
										const std::string& seq1,
										const std::string& seq2,
										const t_loci& prot_loci)
{
	// prot loci to cdna loci
	t_loci loci = prot_loci;
	for_each (beg_to_end(loci), op_prot_loci_to_dna());

	// share processing with protein function
	extract_seq_pair_with_loci(sequences, seq1, seq2, loci);
}

//
//	loci code
//
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	percentage_id

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


//________________________________________________________________________________________

//	get_all_regions_percentage_id

//________________________________________________________________________________________
//	helper functions
struct t_count_identical
{
	explicit t_count_identical(unsigned i = 0): count(i){}
	unsigned count;
	void operator() (char a, char b)
	{
		if (a == b)
			++count;
	}
};

void get_all_regions_percentage_id(vector<string>& proteins, vector<vector<double> >& results)
{
	// percentage id = -1.0 for empty strings
	for (unsigned i = 0; i < proteins.size(); i += 2)
	{
		if (proteins[i].length())
			results[i/2].push_back(
							 for_each_pair(	beg_to_end(proteins[i]),
									proteins[i+1].begin(),
									t_count_identical()).count
									*100.0 / proteins[i].length());
//		cout << "identical =\t" <<
//				for_each_pair(	beg_to_end(proteins[i]),
//									proteins[i+1].begin(),
//									t_count_identical()).count
//			<< "\tout of =\t" << proteins[i].length() << "\t";
	}
}

char table[256][256][256] ;

inline char& table_lookup(char a, char b, char c)
{return table[(unsigned char) a][(unsigned char) b][(unsigned char) c];
}
void fourfoldsymmetry(char table[256][256][256], char nuc1, char nuc2, char value)
{
	table_lookup(nuc1, nuc2, 'A') = value;
	table_lookup(nuc1, nuc2, 'G') = value;
	table_lookup(nuc1, nuc2, 'C') = value;
	table_lookup(nuc1, nuc2, 'T') = value;
}

void twofoldsymmetry_ct(char table[256][256][256], char nuc1, char nuc2, char value)
{
	table_lookup(nuc1, nuc2, 'C') = value;
	table_lookup(nuc1, nuc2, 'T') = value;
}
void twofoldsymmetry_ag(char table[256][256][256], char nuc1, char nuc2, char value)
{
	table_lookup(nuc1, nuc2, 'A') = value;
	table_lookup(nuc1, nuc2, 'G') = value;
}

bool initialize_lookup_table()
{
	for (unsigned  a = 0; a < 256; a++)
		for (unsigned  b = 0; b < 256; b++)
			for (unsigned c = 0; c < 256; c++)
				table [a][b][c] = '\0';
	fourfoldsymmetry(table, 'G', 'C',    1);
	fourfoldsymmetry(table, 'G', 'G',    2);
	fourfoldsymmetry(table, 'C', 'T',    3);	// leucine
	twofoldsymmetry_ag(table, 'T', 'T',  3);	// leucine
	fourfoldsymmetry(table, 'C', 'C',    4);
	fourfoldsymmetry(table, 'C', 'G',    5);	// arginine
	twofoldsymmetry_ag(table, 'A', 'G',  5);	// arginine
	fourfoldsymmetry(table, 'T', 'C',    6);
	fourfoldsymmetry(table, 'A', 'C',    7);
	fourfoldsymmetry(table, 'G', 'T',    8);
	fourfoldsymmetry(table, 'A', 'T',    9);
	twofoldsymmetry_ct(table, 'T', 'G', 10);
	twofoldsymmetry_ct(table, 'G', 'A', 11);
	twofoldsymmetry_ct(table, 'T', 'T', 12);
	twofoldsymmetry_ct(table, 'C', 'A', 13);
	twofoldsymmetry_ct(table, 'A', 'A', 14);
	twofoldsymmetry_ct(table, 'A', 'G', 15);
	twofoldsymmetry_ct(table, 'T', 'A', 16);
	twofoldsymmetry_ag(table, 'G', 'A', 17);
	twofoldsymmetry_ag(table, 'A', 'A', 18);
	twofoldsymmetry_ag(table, 'C', 'A', 19);
	table_lookup('A', 'T', 'G') = 0;		//methionine
	return true;
}
void get_all_regions_cdna_percentage_id(vector<string>& s,
										vector<vector<double> >& ka_results,
										vector<vector<double> >& ks_results)
{
	static bool initialized = initialize_lookup_table();
	// percentage id = -1.0 for empty strings
	for (unsigned i = 0; i < s.size(); i += 2)
	{
		assert(i +1 < s.size());
		if (s[i].length())
		{
			unsigned ka = 0, ks = 0;
			for (unsigned j = 0; j < s[i].length(); j += 3)
			{
				assert(j +2 < s[i].length());
				unsigned cnt_diff = 0;
				if (s[i][j] != s[i+1][j])
					++cnt_diff;
				if (s[i][j+1] != s[i+1][j+1])
					++cnt_diff;
				if (s[i][j+2] != s[i+1][j+2])
					++cnt_diff;
				if (cnt_diff == 1					&&
					table_lookup(s[i][j], s[i][j+1], s[i][j+2]) ==
					table_lookup(s[i+1][j], s[i+1][j+1], s[i+1][j+2]))
					++ks;
				else
					ka += cnt_diff;
			}
			ka_results[i/2].push_back(ka * 100.0 / s[i].length());
			ks_results[i/2].push_back(ks * 100.0  / s[i].length());
		}
	}
	initialized = !initialized;
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	orthologs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888




struct ortholog_region
{
	std::string				gene_id1;
	std::string				gene_id2;
	std::vector<t_segment>	loci_masked;
	std::string				description;
};
void get_gene_ids(std::vector<ortholog_region>& orthologs)
{
	std::string line;
	t_progress_indicator dots(cerr, 10);

	while (cin && getline(cin, line))
	{
		std::vector<std::string> data;
		// get maximum of four parameters
        perl_split (line, "\t", data, false, 4);
		if (data.size() < 4)
		{
			cerr << "\nData not four parameters: " << line << "[" << data.size() << "]\n";
			continue;
		}
		orthologs.push_back(ortholog_region());

		// save gene ids and description
		orthologs.back().gene_id1		= data[0];
		orthologs.back().gene_id2		= data[1];
		orthologs.back().description	= data[3];
		if (! data[2].length())
			continue;
		// split on commas and put into region
		std::vector<std::string> loci_str;
		perl_split(data[2],  "|",  loci_str);
		for (unsigned i = 0; i < loci_str.size(); ++i)
			orthologs.back().loci_masked.push_back(t_segment(loci_str[i]));


		++dots;
	}
}

bool expand_gaps(const string& gaps_str, string& prot, string& cdna, const string& gene_id)
{
	// empty string encodes no gaps!
	if (! gaps_str.length())
		return true;

	// split string of gaps
	vector<string> gaps;
	perl_split(gaps_str,  "|",  gaps);
	for (int i = gaps.size()-1; i >= 0; --i)
	{
		// gaps should not come after a protein or be humungously long
		unsigned beg = atoi(gaps[i].c_str());
		unsigned len = atoi(gaps[i].c_str() + gaps[i].find(':') + 1);
		if (beg > prot.length() || len > 80000)
		{
			return false;
		}
		prot.insert(beg, string(len,  '-'));
		cdna.insert(beg * 3, string(len * 3,  '-'));
	}
	return true;
}

//________________________________________________________________________________________

//	store_gene_ids_as_temp_table

//________________________________________________________________________________________
bool store_gene_ids_as_temp_table(std::vector<ortholog_region>& orthologs, t_pg_cursor& pg_db)
{
	if (!orthologs.size())
		return true;

	cerr << "Give to PostgreSQL the list of gene_ids whose sequence alignments are needed.\n";
	if (!pg_db.ExecCommandOk("CREATE TABLE t_gene_id(gene_id1 TEXT, gene_id2 TEXT, array_index OID)"))
	{
		cerr	<< "\nERROR!! There was an error CREATEing the table of gene_ids "
				<< "t_gene_id\n"
				<<	pg_db.ErrorMessage()
				<<	std::endl;
		return false;
	}

	t_pg_copy_guard copy_buffer(pg_db, "t_gene_id", cerr);
	if (!copy_buffer)
		return false;

	{
		// loop through each sequence
		t_progress_indicator dots(cerr, orthologs.size() / 50 + 1);

		for(unsigned i = 0; i < orthologs.size(); ++i, ++dots)
		{
			copy_buffer.buf += orthologs[i].gene_id1 + "\t" + orthologs[i].gene_id2 + "\t";
			ultoa(i, copy_buffer.buf) += "\n";
			copy_buffer.cached_write();
		}

		if (!copy_buffer.end_copy())
			return false;
	}

	if(!pg_db.ExecCommandOk("CREATE INDEX i_gene_id on t_gene_id(gene_id1)"))
	{
		cerr	<< "\nERROR!! There was an error CREATEing the index for gene_ids "
				<< "t_gene_id\n"
				<<	pg_db.ErrorMessage()
				<<	std::endl;
		return false;
	}

	return true;
}


struct t_error_dump_values
{
	const string& gene_id1;     const string& gene_id2;
	const string& prot1;        const string& prot2;
	const string& cdna1;        const string& cdna2;
	const string& gaps1;        const string& gaps2;
	unsigned align_beg;         unsigned align_end;
	t_error_dump_values(const string& gene_id1_,		const string& gene_id2_,
						const string& prot1_,           const string& prot2_,
						const string& cdna1_,           const string& cdna2_,
						const string& gaps1_,           const string& gaps2_,
						unsigned align_beg_,            unsigned align_end_):
						gene_id1(gene_id1_),            gene_id2(gene_id2_),
						prot1(prot1_),              	prot2(prot2_),
						cdna1(cdna1_),              	cdna2(cdna2_),
						gaps1(gaps1_),                  gaps2(gaps2_),
						align_beg(align_beg_),          align_end(align_end_){}

};
bool dump_params_to_errfile(ofstream& err_strm, const string& err_msg,
							const t_error_dump_values& dump)
{
	//length problem dump diagnostic
	err_strm<< "Error!!!\n\t"
			<< err_msg
			<< "This ortholog will be ignored:\n\t"
			<< "gene_id [ "
							<< dump.gene_id1		<< " / "
							<< dump.gene_id2		<< "]\n\t"
			<< "Gap strings: [ "
							<< dump.gaps1			<< " / "
							<< dump.gaps2			<< "]\n\t"
			<< "Protein lengths: [ "
							<< dump.prot1.length()	<< " / "
							<< dump.prot2.length()	<< "]\n\t"
			<< "Protein 1: [ " << dump.prot1		<< "]\n\t"
			<< "Protein 2: [ " << dump.prot2		<< "]\n\t"
			<< "cDNA lengths: [ "
							<< dump.cdna1.length()	<< " / "
							<< dump.cdna2.length()	<< "]\n\t"
			<< "cDNA 1: [ " << dump.cdna1			<< "]\n\t"
			<< "cDNA 2: [ " << dump.cdna2			<< "]\n\t"
			<< "Alignment [ "
							<< dump.align_beg		<< " / "
							<< dump.align_end		<< "]\n";
	return true;
}


//________________________________________________________________________________________

//	analyze_ortholog

//		die if return false

//________________________________________________________________________________________
bool analyze_ortholog(      std::vector<ortholog_region>& vec_orthologs, ofstream& err_strm,
							vector<vector<double> >& prot_results,
							vector<vector<double> >& cdna_results,
							vector<vector<double> >& ka_results,
							vector<vector<double> >& ks_results,
							unsigned mask_region,
							t_pg_cursor& pg_db, unsigned row)
{
	string prot1 = pg_db.getvalue(row, 0);
	string prot2 = pg_db.getvalue(row, 1);
	string cdna1 = pg_db.getvalue(row, 2);
	string cdna2 = pg_db.getvalue(row, 3);
	string gaps1_str = pg_db.getvalue(row, 4);
	string gaps2_str = pg_db.getvalue(row, 5);
	unsigned align_beg = pg_db.getvalue<unsigned>(row, 6) - 1;
	unsigned align_end = pg_db.getvalue<unsigned>(row, 7) + align_beg;
	unsigned index = pg_db.getvalue<unsigned>(row, 8);

	if (index > vec_orthologs.size())
	{
		cerr	<< "\n\nError: The index returned [" << index << "]"
				<< " was greater than the rows of gene_ids [" << vec_orthologs.size() <<"]\n";
		return false;
	}
	ortholog_region& orthologs = vec_orthologs[index];

	// check gaps, alignment and lengths
	t_error_dump_values dump (orthologs.gene_id1,   orthologs.gene_id2,
							  prot1,                prot2,
							  cdna1,                cdna2,
							  gaps1_str,            gaps2_str,
							  align_beg,            align_end);
	if (!expand_gaps(gaps1_str, prot1, cdna1, orthologs.gene_id1) ||
		!expand_gaps(gaps2_str, prot2, cdna2, orthologs.gene_id2))
		return dump_params_to_errfile(err_strm, "Errors in Gap Lengths. ", dump);
	if (align_beg >= align_end || align_end > prot1.length())
		return dump_params_to_errfile(err_strm, "Alignment coordinates invalid. ", dump);
	if (prot1.length() != prot2.length() ||
		cdna1.length() != cdna2.length() ||
		cdna1.length() != prot1.length() * 3)
		return dump_params_to_errfile(err_strm, "Unmatched sequence strings. ", dump);

	//
	// merge overlapping regions
	//
	std::vector<t_loci> vec_loci(3);
	vec_loci[0].push_back(t_segment(align_beg, align_end));
	vec_loci[1].swap(orthologs.loci_masked);
	merge_overlapping_loci(vec_loci[1], vec_loci[2], align_beg, align_end);

	if (!(mask_region & eINVERSE))
		vec_loci.erase(vec_loci.begin() + 2);
	if (!(mask_region & eMASKED))
		vec_loci.erase(vec_loci.begin() + 1);
	if (!(mask_region & eWHOLE))
		vec_loci.erase(vec_loci.begin() + 0);
	assert(vec_loci.size());


	//____________________________________________________________________________________
	// get strings comprising concatenated domain/free regions
	//____________________________________________________________________________________
		vector<string> proteins;
		for (unsigned j = 0; j < vec_loci.size(); ++j)
			extract_seq_pair_with_loci(proteins, prot1, prot2, vec_loci[j]);

	//____________________________________________________________________________________
	// get protein percentage ids
	//____________________________________________________________________________________
		get_all_regions_percentage_id(proteins, prot_results);



	//____________________________________________________________________________________
	// get nucleotide percentage ids
	//____________________________________________________________________________________
		vector<string> cDNAs;
		for (unsigned j = 0; j < vec_loci.size(); ++j)
			extract_cdna_pair_with_loci(cDNAs, cdna1, cdna2, vec_loci[j]);
		get_all_regions_percentage_id(cDNAs, cdna_results);
		get_all_regions_cdna_percentage_id(cDNAs, ka_results, ks_results);

//		cout << "\n";

	return true;
}

void print_median(std::vector<double>& results)
{
	if (!results.size())
		return;
	std::sort(beg_to_end(results));
	unsigned half		= int(results.size() * 0.5);
	unsigned sixteen	= int(results.size() * 0.166666666666666);
	unsigned eightythree= int(results.size() * 0.833333333333333);
	cout 	<<	std::setprecision(3)
		 	<<	std::fixed
			<< "\tmedian =\t"
			<< results[half]	<< "%\t[ "
			<< results[sixteen] << "% - "
			<< results[eightythree] << "% ]\n";
}

//________________________________________________________________________________________

//	analyse_all_orthologs

//		die if return false

//________________________________________________________________________________________
bool analyse_all_orthologs(std::vector<ortholog_region>& orthologs, const string& sql_str,
							unsigned cnt_estimate,
							unsigned mask_region,
							const std::string& error_file,
							t_pg_cursor& pg_db)
{
	ofstream err_strm (error_file.c_str());
	if (!err_strm)
	{
		cerr << "\n\nError\n\tUnable to open " << error_file << "\n";
		return false;
	}


	// create binary cursor
	if (!pg_db.Declare(sql_str.c_str(), true))
	{
		cerr << "Failed to create a cursor to retrieve orthologuous sequences!\n";
		return false;
	}


	vector<vector<double> > prot_results(3);
	vector<vector<double> > cdna_results(3);
	vector<vector<double> > ka_results(3);
	vector<vector<double> > ks_results(3);

	if (mask_region & eINVERSE)
	{
		cdna_results.pop_back();
		prot_results.pop_back();
		ka_results.pop_back();
		ks_results.pop_back();
	}
	if (!(mask_region & eMASKED))
	{
		cdna_results.pop_back();
		prot_results.pop_back();
		ka_results.pop_back();
		ks_results.pop_back();
	}
	if (!(mask_region & eWHOLE))
	{
		cdna_results.pop_back();
		prot_results.pop_back();
		ka_results.pop_back();
		ks_results.pop_back();
	}


	// time whole process
	t_portable_timer timer;
	{
		// get data in batches
		t_progress_indicator dots(cerr, 1);

		unsigned batch_size = (cnt_estimate / 50 + 1);
		while (pg_db.Fetch(batch_size))
		{
			int curr_batch_size = pg_db.Tuples();
			if (!curr_batch_size)
				break;
			dots++;

			// buffer output
			for (int i = 0; i < curr_batch_size; ++i)
				if (!analyze_ortholog(orthologs, err_strm,
										prot_results, cdna_results,
										ka_results, ks_results,
										mask_region,
										pg_db, i))
					return false;
		}
	}
	double seconds = timer.reset();
	if (seconds > 120)
		cerr << "\tAnalysing orthologs took\t" << seconds / 60 << " minutes\n";
	else
		cerr << "\tAnalysing orthologs took\t" << seconds << " seconds\n";

	// close cursor
	if (!pg_db.Close())
		return pg_print_sql_err("close a cursor for retrieving ortholog sequences",
									pg_db, "CLOSE panda_cursor", cerr);

	unsigned results_pos = 0;
	if (mask_region & eWHOLE)
	{
		cout << "Whole gene %Protein Id\n";
		print_median(prot_results[results_pos]);
		cout << "Whole gene %cDNA Id\n";
		print_median(cdna_results[results_pos]);
		cout << "Synonymous changes %cDNA Id\n";
		print_median(ks_results[results_pos]);
		cout << "Non-synonymous changes %cDNA Id\n";
		print_median(ka_results[results_pos]);
		++results_pos;
	}
	if (mask_region & eMASKED)
	{
		cout << "Masked %Protein Id\n";
		print_median(prot_results[results_pos]);
		cout << "Masked %cDNA Id\n";
		print_median(cdna_results[results_pos]);
		cout << "Synonymous changes %cDNA Id\n";
		print_median(ks_results[results_pos]);
		cout << "Non-synonymous changes %cDNA Id\n";
		print_median(ka_results[results_pos]);
		++results_pos;
	}
	if (mask_region & eINVERSE)
	{
		cout << "Unmasked %Protein Id\n";
		print_median(prot_results[results_pos]);
		cout << "Unmasked %cDNA Id\n";
		print_median(cdna_results[results_pos]);
		cout << "Synonymous changes %cDNA Id\n";
		print_median(ks_results[results_pos]);
		cout << "Non-synonymous changes %cDNA Id\n";
		print_median(ka_results[results_pos]);
		++results_pos;
	}
	return true;
}










//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	// don't use C streams
	std::ios::sync_with_stdio(false);


	// get command line arguments

	// if no arguments print usage and exit
	if (argc == 1)
	{
		char * my_argv[2] ={ argv[0], "-h"};
		int argc = 2;
		t_program_args args;
		if (!process_command_line_options(argc, my_argv, args))
			return 1;
	}
	t_program_args args;
	if (!process_command_line_options(argc, argv, args))
		return 1;


	//____________________________________________________________________________________
	// 		Connect to PostgreSQL
	//____________________________________________________________________________________
	t_pg_cursor pg_db(args.db_login, "panda_cursor");	//"dbname=test user=postgres"
	if ( pg_db.ConnectionBad() )
	{
		cerr		<< "\nERROR!! Connection to database ["
					<< args.db_login << "] failed.\n"
					<< pg_db.ErrorMessage() << "\n\n";
		return 1;
	}



	//____________________________________________________________________________________
	// read ortholog data from STDIN
	//____________________________________________________________________________________
	std::vector<ortholog_region> orthologs;
	std::cerr << "Read supplied gene IDs...\n";
	get_gene_ids(orthologs);

	std::cerr << "Store as temp table...\n";
	store_gene_ids_as_temp_table(orthologs, pg_db);

	std::cerr << "Count number of orthologs sequences to be retrieved...\n";
	//____________________________________________________________________________________
	// analyse each ortholog, retrieving the sequence from the specified database
	//____________________________________________________________________________________
	string sql_fields = "SELECT s1.sequence, s2.sequence, cdna1.cdna, cdna2.cdna, "
							"aln_gaps1, aln_gaps2, aln.bl2seq_align_start, "
							"aln.bl2seq_align_len, array_index ";
	string sql_tables = "FROM t_gene_id NATURAL JOIN orthologs." + args.taxon1 + '_' +
						args.taxon2 + "_alignment aln JOIN " +
						args.taxon1 + ".cdna_sequences cdna1 ON (aln.pi1 = cdna1.pi AND "
						"aln.acc_code1 = cdna1.acc_code) JOIN sequences s1 ON "
						"(aln.pi1 = s1.pi) JOIN " + args.taxon2 + ".cdna_sequences cdna2 "
						"ON (aln.pi2 = cdna2.pi AND aln.acc_code2 = cdna2.acc_code) JOIN "
						"sequences s2 ON (aln.pi2 = s2.pi)";

	if (!pg_db.ExecTuplesOk(("SELECT count (*) " + sql_tables).c_str()))
	{
		cerr << "Error:\n\tCould not retrieve the count of ortholog sequences to be "
					"aligned.\n\n";
		return 1;
	}
	unsigned cnt_estimate =	atoi(pg_db.getvalue(0, 0));
	cerr << "\t" << cnt_estimate << "\tortholog sequences to be processed...\n";



	if (!analyse_all_orthologs(orthologs, sql_fields + sql_tables, cnt_estimate, args.region,
								args.error_file,
								pg_db))
	{
		cerr << "Analysis failed ;-<\n";
		return 1;
	}


	cerr << "Success!\n";
	return 0;
}


